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ABSTRACT 


The  three-dimensional  shear-driven  cavity  flow  is  numerically  investigated  at  Reynolds 
numbers  of  5000  and  10,000.  This  investigation  focuses  on  the  unsteadiness  and  turbulent  char¬ 
acteristics  of  the  flow.  At  the  moderate  Reynods  number  (Re  =  5000)  where  the  cavity  flow  is 
fully  laminar,  a  direct  numerical  simulation  (DNS)  is  used;  at  the  higher  Reynolds  number 
(Re  =10,000),  large-eddy  simulation  (LES)  methodology  is  used  to  predict  the  cavity  flow. 
Establishing  a  suitable  form  for  the  subgrid  scale  (SGS)  turbulence  model  in  this  complex  flow  is 
guided  by  the  DNS  results  at  Re  =  5000.  Additionally,  the  SGS  model  is  verified  against  DNS 
results  at  Re  =  7500,  where  the  cavity  flow  is  known  through  experimentation  to  be  locally 
transitional.  The  LES  results  verify  the  published  experimental  evidence,  as  well  as  introduce 
new  flow  features  within  the  cavity. 
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INTRODUCTION 


For  more  than  three  decades,  the  shear-driven  cavity  flow  problem  has  served  as  an  excel¬ 
lent  test  case  for  verifying  new  or  improved  numerical  solution  techniques  of  incompressible 
flows.  Usually,  the  test  case  involves  simulating  a  two-dimensional  (2D)  cavity  flow  at  a  low 
Reynolds  number  (Re  <  1000).  Under  these  conditions,  the  flow  is  strictly  laminar  and  steady. 
The  primary  purpose  of  the  simulations  is  to  illustrate  the  rapid  speed  of  convergence  to  steady 
state  and  the  solution  method’s  ability  to  capture  the  basic  features  of  the  flow.  Some  2D  simu¬ 
lations  (for  example,  Ghia  et  al.  [1]  and  Gustafson  and  Halasi  [2])  include  a  discussion  of  the 
flow  characteristics  and  have  revealed  important  salient  features  of  the  steady  flow  at  much 
higher  Reynolds  numbers  (Re  <  10,000).  The  extensive  results  of  Ghia  et  al.,  in  particular,  serve 
most  often  as  a  base  for  comparison  of  new  2D  predictions  because  of  the  fine  grid  resolution 
they  used  at  the  corresponding  Re.  Only  a  few  numerical  studies  have  reported  the  turbulent 
characteristics  of  the  shear-driven  cavity  flow  (Gosman  et  al.  [3],  Young  et  al.  [4],  Ideriah  [5], 
and  Gaskell  and  Lau  [6],  for  instance).  Each  study  simulated  a  2D  geometry  with  the  turbulence 
fully  modeled.  Treatment  of  the  results  was  qualitative,  meaning  that  the  investigators  validated 
the  particular  numerical  scheme  or  evaluated  various  turbulence  models.  Separate  studies  show¬ 
ing  the  flow  evolution  in  the  2D  cavity  under  an  impulsively  started  lid  and  an  oscillating  lid 
were  reported  by  Soh  and  Goodrich  [7].  In  both  simulations,  the  flow  was  laminar  (Re  =  400) 
and  the  predictions  continued  until  steady-state.  Published  results  from  three-dimensional  (3D) 
simulations  include  those  of  Moin  and  Kim  [8],  Freitas  et  al.  [9],  and  Prasad  et  al.  [10],  where 
the  Reynolds  number  was  restricted  to  low-to-moderate  values  (Re  <  3200).  There,  the  flow  was 
also  laminar.  The  simulations  showed  the  appearance  of  pairs  of  quasi-steady  and  unsteady 
spanwise  Taylor-Gortler-like  (TGL)  vortices  along  the  cavity  bottom.  Besides  validating  the 
conservativeness  of  their  particular  solution  technique,  each  revealed  this  important  3D  charac¬ 
teristic,  which  had  been  observed  experimentally  but  not  numerically. 

In  addition  to  furnishing  us  with  a  classic  problem  for  validating  solution  techniques  of 
incompressible  flows,  the  predominant  features  of  the  shear-driven  cavity  flow  also  have  impor¬ 
tant  physical  significance  in  view  of  engineering  design.  For  example,  the  main  flow  character¬ 
istics  created  by  3D  geometries  such  as  recessed  cavities,  recurrent  ribs  or  slots  and  curved  ducts, 
for  the  purposes  of  mass  and/or  energy  exchange,  are  similar  to  those  of  this  model  problem. 
Given  this  fact,  the  objective  of  the  present  work  is  to  identify  and  characterize,  through  numeri¬ 
cal  simulation,  the  unsteadiness  and  turbulence  in  a  3D  shear-driven  cavity  flow.  This  numerical 
investigation  is  presented  for  Reynolds  numbers  of  5000  and  10,000.  Numerical  results  of  the 
unsteady  and  laminar  flow  characteristics  at  much  lower  Reynolds  numbers 
(Re  =  2000  and  Re  =  3200)  were  presented  by  Jordan  and  Ragab  [11]. 

For  the  moderate  Reynolds  number  test  case  (Re  =  5000)  where  the  flow  is  unsteady  and 
laminar,  the  computation  was  a  direct  numerical  simulation  (DNS).  At  the  higher  Re,  the  predic¬ 
tion  method  was  a  large-eddy  simulation  (LES).  Although  extensive  experimental  results  of  the 
high-Re  cavity  flow  have  been  reported  for  just  over  10  years  now,  this  LES  investigation  is  a 
first  attempt  to  study  the  turbulent  flow  characteristics  numerically.  Before  deriving  any  new 
unsteady  or  turbulent  flow  features,  however,  the  numerical  results  are  first  scrutinized  thor- 
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oughly  against  the  experimental  evidence.  Thus,  this  LES  computation  verifies  the  published 
experimental  data  as  well  as  introduces  new  characteristics  about  the  flow. 


BASIC  MODEL  PROBLEM  INCLUDING  EXPERIMENTAL  OBSERVATIONS 

Geometrically,  the  classic  2D  problem  consists  of  a  closed  cavity  of  unit  height  (H)  and 
unit  width  (W)  with  a  lid  moving  horizontally  at  unit  velocity  (U).  The  Reynolds  number  is 
defined  as  Re  =  UW/v,  where  v  is  the  kinematic  viscosity.  The  associated  recirculation  flow  is 
characterized  basically  by  a  primary  vortex,  a  downstream  secondary  eddy,  an  upstream  second¬ 
ary  eddy,  and  an  upper  secondary  eddy.  These  basic  features  are  sketched  in  figure  1.  The  upper 


Figure  1.  Sketch  of  the  Basic  Features  of  Recirculation  in  the  Two-Dimensional 

Shear-Driven  Cavity  Flow  Problem 


secondary  eddy  appears  at  Re  >  3200.  At  higher  Re  (>  5000)  the  2D  predictions  show  a  tertiary 
eddy  in  each  of  the  lower  corners.  Based  on  the  experimental  observations  of  Koseff  and  Street 
[12]  of  the  3D  cavity  flow,  these  tertiary  eddies  are  unique  only  to  the  2D  simulations.  From 
their  observations,  Koseff  and  Street  also  reported  that  the  basic  features  of  the  recirculation  flow 
in  figure  1  are  nearly  self-similar  throughout  the  cavity  span  at  high  Reynolds  numbers  where  the 
flow  is  locally  transitional. 

In  the  3D  cavity  new  vortical  structures  are  formed.  In  a  spanwise  plane,  there  are  several 
pairs  of  TGL  vortices  and  a  lower  corner  vortex  at  the  end-walls  (see  sketch  in  figure  2).  Ac¬ 
cording  to  Koseff  and  Street  [13]  and  Prasad  et  al.  [10],  the  impetus  for  manifesting  the  TGL 
vortices  is  the  instability  of  the  concave  free  shear  layer  that  separates  the  primary  vortex  from 
the  downstream  secondary  eddy.  Generation  of  the  vortices  occurs  just  above  the  concave 
surface,  much  like  the  experimental  observations  of  Taylor  [14]  for  the  flow  between  rotating 
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Figure  2.  Sketch  of  the  Taylor-Gortler-Like  Vortex  Pairs  as  Observed  in  the  Flow 

Visualization  Experiments  [12] 


cylinders  and  also  the  concave  boundary  layer  investigated  by  Gortler  [15].  The  size  and  number 
of  pairs  of  TGL  vortices  depends  strongly  on  the  Reynolds  number  and  the  cavity  spanwise 
aspect  ratio  (SAR).  In  the  flow  visualization  results  reported  by  Rhee  et  al.  [16]  for  Re  <  6000 
and  SAR  =  3.0,  the  spanwise  flow  maintained  symmetry  about  the  mid-span  plane.  Although  the 
TGL  vortices  meander  slowly  along  the  cavity  bottom  at  moderate  Re,  Koseff  and  Street  [1 3] 
noted  that  their  basic  spanwise  flow  character  still  remains  symmetric.  The  other  important 
feature  in  the  spanwise  direction  is  the  lower  corner  vortex.  The  origin  of  this  flow  structure  was 
explained  by  Koseff  and  Street  [17]  after  they  examined  the  experimental  results  of  de 
Brederoede  and  Bradshaw  [18].  Manifestation  of  this  vortex  is  a  consequence  of  the  shear  and 
pressure  force  adjustment  in  the  streamwise  recirculating  flow  caused  by  the  no-slip  condition 
along  the  spanwise  end-wall.  Like  the  TGL  vortices,  the  comer  vortex  becomes  unsteady  at 
moderate  Reynolds  numbers  (Re  >  3200).  The  size  and  extent  of  the  comer  vortex  strongly 
influences  the  TGL  vortices.  Thus,  the  numerical  simulations  must  provide  sufficient  resolution, 
spatially  and  temporally,  to  capture  its  characteristics  accurately. 

Experimental  observations  show  the  first  sign  of  turbulence  taking  place  within  the  free 
shear  layer  that  lies  between  the  primary  vortex  and  downstream  secondary  eddy.  This  local 
transition  to  turbulence  occurs  at  a  Reynolds  number  somewhere  between  6000  and  8000  [13]. 
The  flow  within  that  region  is  unsteady.  If  the  Reynolds  number  is  increased,  the  onset  of 
turbulence  diffuses  the  TGL  vortices  such  that  their  deterministic  structure  becomes  obscured. 

At  Re  =  10,000,  frequency  spectra  of  both  the  horizontal  and  vertical  fluctuations  within  the 
region  of  the  free  shear  layer  display  an  inertial  subrange.  It  should  be  noted  that  the  visualiza¬ 
tion  results  show  the  flow  within  the  free  shear  layer  at  Re  =  10,000  as  being  still  transitional. 
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GOVERNING  EQUATIONS  AND  METHOD  OF  SOLUTION 


The  governing  LES  equations  of  the  resolvable  field  are  the  spatially  filtered  incompress¬ 
ible  Navier-Stokes  equations: 


du,  d  / _ \  op 


dp  1  d2u.  dtjj 


dxj 


dx, 


+ 


Re  dxjdxj 


dxj 


(1) 


3x, 


(2) 


where  p  is  the  pressure,  and  tjj  is  the  subgrid  scale  stress  tensor,  which  is  defined  as 


S  =  u,Uj  -  u,Uj . 


The  system  of  equations  is  time-advanced  according  to  a  variant  of  the  fractional-step 
method  [19].  The  Crank-Nicolson  scheme  is  applied  to  the  viscous  terms  to  eliminate  the  vis¬ 
cous  stability  restriction.  An  explicit  three-step  third-order-accurate  Runge-Kutta  (R-K)  scheme 
is  used  for  the  convective  terms  instead  of  the  more  popular  Adams-Bashforth  method  (see  Clark 
et  al.  [20],  for  example).  This  is  done  because  the  Adams-Bashforth  method  is  weakly  unstable 
when  applied  to  the  linear  convection  equation.  Also,  the  explicit  form  of  the  R-K  scheme  easily 
permits  high-order  spatial  differences.  In  this  application  of  the  LES  approach,  the  convective 
terms  are  spatially  discretized  by  third-order  upwind-biased  finite  differences,  while  the  diffusive 
terms  are  differenced  by  second-order  central  differences.  At  the  points  next  to  the  wall  bound¬ 
aries,  exterior  points  necessary  to  complete  the  upwind  differencing  are  obtained  through  ex¬ 
trapolation  of  the  interior  field  points.  All  terms  in  Poisson’s  equation  for  solution  of  the  pres¬ 
sure  variable  are  central  differenced  to  the  second  order.  Strong  coupling  between  the  pressure 
and  velocity  components  is  maintained  through  a  fourth-order-accurate  compact  differencing 
scheme  for  the  pressure  gradient  in  the  velocity  update  equation.  This  discretization  scheme 
provides  an  overall  method  that  is  second-order-accurate  in  both  space  and  time. 

The  velocity  components  are  collocated  with  the  grid  points  to  permit  easy  development  of 
a  set  of  wall  boundary  conditions  for  the  velocity  field.  Approximate-factorization  (A-F)  is 
applied  for  solution  of  the  intermediate  velocity  components.  The  pressure  field  is  staggered 
from  the  velocity  field  to  eliminate  spurious  oscillations  in  the  flow  solutions  [21],  No  boundary 
condition  is  required  for  the  pressure  variable.  Poisson’s  equation  for  solution  of  the  pressure 
variable  is  recast  into  a  residual  form  and  solved  by  a  variant  of  the  modified  strongly  implicit 
procedure  [22].  This  procedure  is  implemented  as  an  implicit  elliptic  solver  of  a  planar  surface 
that  moves  through  the  volume.  After  the  initial  time  steps,  only  a  few  iterations  on  the  pressure 
variable  are  necessary  to  satisfy  the  incompressibility  constraint  at  each  of  the  following  time 
steps.  Extensive  details  of  these  numerics,  along  with  several  test  cases,  can  be  found  in  Jordan 
and  Ragab  [11]. 
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SUBGRID  TURBULENCE  MODEL 


For  the  computations  of  the  turbulent  shear-driven  cavity  flow,  the  Smagorinsky  eddy 
viscosity  model  [23]  was  implemented  to  represent  the  subgrid  scale  (SGS)  turbulence.  In  tensor 
notation  the  Smagorinsky  model  is: 


^ij  3  ^ij^kk  —  2vTSy,  (3) 

Sij  =i(uiJ+uJ>1),  (4) 


VT  —  L  ^SySjj  , 


(5) 


L  =  ct 


l-exp(-y+/A+)m|  (AjA2A3) 


(6) 


In  this  model,  L  is  the  turbulent  characteristic  length  scale.  Modification  of  the  length  scale  by  a 
form  of  van  Driest  damping  [24]  was  necessary  to  account  for  the  effects  of  the  solid  wall 

boundaries.  The  length  parameter  y+  is  the  minimum  field  value  of  y+  =  y^/x,  /  P  /  v,  where 

is  the  magnitude  of  the  local  shear  stress.  The  subscripts  1,  2,  and  3  of  the  filter  width  D 

permit  versatile  filtering  for  modeling  anisotropic  flows.  As  a  minimum  A,  =  2h1;  where  h,  is 
the  grid  spacing  in  the  i  direction.  The  exponents  m  and  n  in  the  van  Driest  damping  function 
were  determined  by  analyzing  the  DNS  results  and  by  qualitative  comparisons  to  the  published 
experimental  measurements.  The  importance  of  including  van  Driest  damping  in  the 
Smagorinsky  model  became  apparent  after  applying  the  model  to  the  DNS  results  at  Re  =  5000 
without  the  damping  modification.  The  model  repeatedly  showed  high  levels  of  turbulent  eddy 
viscosity  near  the  upper  half  of  the  downstream  wall.  At  times,  these  levels  reached  the  same 
order  of  magnitude  as  the  kinematic  viscosity,  marking  the  onset  of  turbulence.  Based  on  the 
experimental  observations,  however,  this  is  a  false  indication.  As  discussed  earlier,  the  flow 
visualization  results  suggest  that  initial  signs  of  turbulence  occur  within  the  free  shear  layer 
between  the  primary  vortex  and  the  downstream  secondary  eddy  located  near  the  lower  half  of 
the  downstream  wall.  By  damping  the  length  scale  according  to  a  form  of  van  Driest,  the  near¬ 
wall  turbulent  eddy  viscosity  was  reduced  and  the  modified  model  agreed  qualitatively  with  the 
experimental  evidence.  This  issue  is  addressed  further  in  the  next  section. 


RESULTS  AND  DISCUSSION 

This  section  presents  the  unsteady  and  turbulent  flow  results  from  a  numerical  investigation 
of  the  3D  shear-driven  cavity  and  includes  comparisons  with  the  published  experimental  obser¬ 
vations  and  measurements.  This  investigation  is  of  a  cavity  with  SAR  =  3.0.  Reynolds  numbers 
of  5000  and  10,000  were  simulated.  The  cavity  geometry  was  modeled  as  unit  streamwise 
(x-direction)  and  normal  (y-direction)  lengths;  W  =  1.0  and  H  =  1.0  (see  figure  3).  Based  on  the 
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Figure  3.  Model  Problem  for  Present  DNS  and  LES  Predictions  of  the 

Shear-Driven  Cavity  Flow 


published  flow  visualization  data  [13,16,17]  for  Re  <6000  and  previous  simulations  [9,10,1 1]  for 
Re  =  3200,  only  one-half  of  the  cavity  span  was  modeled.  Thus,  the  spanwise  (z-direction) 
length  is  1 .5  times  that  of  the  streamwise  length,  with  one  boundary  modeled  as  a  plane  of 
symmetry  at  mid-span  (z  =  0.0)  and  the  other  a  solid  end-wall  (z  =  1.5).  For  the  Re  =  10,000 
simulation,  a  localized  effect  of  the  symmetry  plane  assumption  on  the  numerical  results  will  be 
shown.  The  lid  moved  horizontally  with  unit  velocity  (Ul=1.0);  hence,  the  Reynolds  number  is 
simply  the  inverse  of  the  kinematic  viscosity.  The  no-slip  condition  was  enforced  along  all 
boundaries  except  at  the  mid-span  plane,  which  was  treated  numerically  as  a  plane  of  symmetry. 
All  grids  selected  had  uniform  point  spacing. 

Each  simulation  was  initialized  by  an  impulsively  started  lid.  Inasmuch  as  the  convective 
terms  are  time-split  by  the  R-K  technique,  extremely  small  time  increments  were  not  necessary 
early  in  the  simulation  to  maintain  stability.  For  both  simulations,  a  low  CFL  value  of  0.5  was 
chosen  to  ensure  proper  temporal  resolution  of  the  flow  unsteadiness  rather  than  controlling 
numerical  stability.  It  was  found  that  CFL  values  near  the  stability  limit  distorted  the  TGL 
vortex  structures.  Herein,  computational  results,  which  are  labeled  t  =  0.0,  represent  solutions 
that  had  been  time-advanced  until  the  effects  of  the  impulsively  started  lid  on  the  flow  evolution 
became  negligible.  In  both  simulations,  this  time  was  approximately  1  minute. 

In  view  of  the  experimental  evidence,  the  cavity  flow  is  entirely  laminar  at  low  to  moderate 
Reynolds  numbers  (Re  <  5000).  Therefore,  the  simulation  at  Re  =  5000  is  essentially  a  DNS 
prediction.  The  importance  of  the  DNS  results  with  respect  to  the  LES  methodology  lies  in 
improving  the  estimate  of  the  turbulence  length  scale;  specifically,  determining  suitable  values 
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for  the  exponents  (m  and  n)  in  the  van  Driest  damping  function  and  for  the  Smagorinsky  constant 
(Cs).  The  procedure  used  here  closely  follows  the  SGS  turbulence  model  development  by  Clark 
et  al.  [20].  The  model  parameters  were  judged  according  to  the  best  overall  set  of  coefficients 
attained  from  correlations  between  the  exact  (computed  from  the  DNS  data)  and  the  model 
results  in  terms  of  the  SGS  stresses  tjj.  The  correlation  coefficient  is  defined  as 

<  Cjj  >=<  ejjiiijj  >  /  <  e,j  >1/2<  m2  >1/2  where  <  ey  >  and  <  >  are  the  spatially  en¬ 

semble-averaged  exact  and  model  tjj  components,  respectively.  Inasmuch  as  the  spanwise  flow 
direction  was  assumed  to  be  homogeneous,  only  the  coefficients  cj  j,  C22>  and  c|2  were  deter¬ 
mined.  In  addition,  only  positive  values  of  these  coefficients  were  considered  in  the  model 
evaluation.  Before  the  model  was  implemented  in  the  Re  =  10,000  simulation,  its  quality  was 
checked  against  the  published  experimental  observations.  In  this  section,  we  will  show  that  the 
definition  in  equation  (6)  for  the  turbulence  length  scale,  along  with  its  associated  constants,  is 
acceptable  for  the  shear-driven  cavity  flow  problem  at  the  simulated  Reynolds  numbers.  Al¬ 
though  designing  a  proper  form  for  the  turbulence  model  was  focused  on  the  DNS  results  at 
Re  =  5000,  the  model  quality  was  also  checked  at  a  higher  Reynolds  number  where  the  flow  is 
locally  transitional. 


DNS  RESULTS 

For  the  test  case  of  Re  =  5000,  DNS  computations  were  performed  over  a  65  x  65  x  65 
(x,  y,  z-direction)  uniform  grid.  Based  on  the  grids  used  in  previous  simulations  of  3D  cavity 
flow  [9,10,1 1],  this  mesh  resolution  is  adequate  for  this  Reynolds  number.  The  simulation  was 
carried  out  for  approximately  4  minutes.  Figures  4a  and  4b  show  a  set  of  snapshots  of  the  un¬ 
steady  flow  results  at  time  t  =  15.0.  This  relative  time  corresponds  to  an  actual  time  of  over  1 
minute  since  impulsive  start  of  the  lid.  The  velocity  vectors  in  figure  4a  represent  the  recircula¬ 
tion  flow  at  the  mid-span  plane,  whereas  those  in  figure  4b  depict  the  spanwise  flow  at  plane  x  = 
0.77.  For  clarity,  each  of  the  velocity  vectors  is  normalized  by  its  own  magnitude.  Although 
there  seems  to  be  3D  effects  in  the  upper  half  of  the  spanwise  plane,  their  magnitudes  are  small 
and  therefore  have  little  significance  for  the  flow  structure.  At  this  instant  in  time  (t  =  15.0),  the 
three  basic  features  of  the  recirculation  flow  that  are  common  to  the  2D  simulations  are  distinctly 
visible.  Likewise,  the  primary  vortex  core  is  positioned  close  to  the  cavity  center.  This  agree¬ 
ment  is  due  primarily  to  the  minimal  influence  of  the  spanwise  flow  on  the  mid-span  recircula¬ 
tion  flow  (see  figure  4b).  The  spanwise  flow  vectors  show  four  TGL  vortices  of  nearly  the  same 
height  that  lie  fully  within  the  cavity  span.  Since  the  static  pressure  attains  a  minimum  within 
the  vortex  core,  contours  of  the  pressure  variable  can  verify  the  existence  of  each  vortex.  An 
example  of  this  is  portrayed  in  figure  4c.  A  few  grid  lines  are  superimposed  over  the  pressure 
contours  and  the  corresponding  velocity  vectors  to  help  identify  each  respective  vortex. 

At  time  t  =  181.0,  the  snapshots  paint  a  very  different  picture  (see  figures  5a  and  5b).  The 
3D  effects  on  the  basic  recirculating  flow  features  are  clearly  displayed.  For  example,  the  TGL 
vortex  that  straddles  the  mid-span  plane  severely  distorts  the  basic  structure  of  the  downstream 
secondary  eddy  (DSE).  Throughout  most  of  the  simulation,  a  TGL  vortex  structure  straddled  the 
mid-span  plane,  which  precluded  development  of  the  local  DSE.  According  to  the  experimental 
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Figure  4.  Snapshots  of  the  Unsteady  Shear-Driven  Cavity  Flow  at  Relative  Time 
t  =  15.0;  (a)  Recirculation  at  Mid-Span,  (b)  TGL  Vortices  at  Plane  x  =  0.77 and 

(c)  Spanwise  Static  Pressure  Contours 


results  [16],  eight  pairs  of  TGL  vortices  were  visualized  at  Re  =  3200,  and  1 1  vortex  pairs  at  Re 
=  6000.  At  Re  =  5000  we  found  nine  vortex  pairs;  one  typically  straddling  the  mid-span  plane 
and  four  others  spanning  the  cavity  floor.  This  result  is  shown  in  figure  6,  which  is  a  plot  of  the 
x-vorticity  at  completion  of  the  simulation  (t  =181.0).  Notice  that  the  streamwise  extent  of  each 
vortex  pair  does  not  strongly  interact  with  the  upstream  secondary  eddy  (USE).  This  is  because 
the  flow  process  of  fluid  entrainment  from  the  primary  vortex  to  sustain  the  structural  integrity  of 
each  TGL  vortex  is  terminated  once  the  primary  vortex  separates  upstream.  Particle  traces  that 
illustrate  this  flow  process  are  presented  in  the  next  section. 

In  the  literature,  neither  velocity  time  traces  nor  mean  velocity  experimental  data  appear  for 
the  3D  cavity  flow  at  Re  =  5000  and  SAR  =  3,0.  Here,  shown  in  figures  7a  through  7d  are  time 
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Figure  6.  Spanwise  Distribution  and  Streamwise  Extent  of  TGL  Vortices  at 
Completion  of  DNS  Computation  (t  =  3  min) 
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traces  of  a  vertical  and  a  horizontal  velocity  component  and  their  power  spectra  for  the 
Re  =  5000  test  case.  The  traces  were  extracted  from  recordings  taken  in  the  vicinity  of  the 
downstream  free  shear  layer.  For  reference,  the  power  spectra  include  Kolmogorov’s  slope  of 
the  inertial  sub-range.  Both  spectra  show  numerous  amplified  frequencies,  signifying  an  un¬ 
steady  flow  that  is  still  deterministic.  In  figure  8  the  computed  mean  horizontal  velocity  along 
the  mid-span  centerline  is  compared  with  the  reported  experimental  data  for  SAR  =1:1  and  SAR 
=  0.5:1  [26],  The  DNS  profile  illustrates  a  further  weakening  of  the  “energy-sink”  effect  of  the 
spanwise  end-walls  on  the  primary  recirculation  vortex  core  when  the  cavity  SAR  is  extended  to 
three. 


As  noted  earlier,  the  DNS  results  of  this  test  case  helped  to  establish  an  acceptable  form 
of  the  SGS  turbulence  model  for  the  LES  predictions.  In  particular,  parametric  studies  were 
performed  on  the  DNS  results  at  intermittent  time  intervals  to  find  a  global  set  of  values  for 
Smagorinsky’s  constant  (Cs)  and  the  exponents  (m  and  n)  in  the  van  Driest  damping  function. 
Furthermore,  signs  of  transition  first  emerge  within  the  downstream  free  shear  layer  [13].  By 
knowing  this  information  a  priori,  it  is  also  possible  to  perform  quality  checks  on  the  model  in 
terms  of  the  turbulent  eddy  viscosity  (TEV)  magnitudes  and  distribution.  In  the  following 
figures,  the  TEV  magnitudes  are  normalized  by  the  kinematic  viscosity.  Since  the  TGL  vortices 
severely  disturbed  the  basic  features  of  the  recirculating  flow,  the  model  quality  was  inspected 
only  on  planes  lying  between  the  spanwise  vortex  pairs. 

The  best  set  of  fully-averaged  correlation  coefficients  attained  from  the  parametric  study 
was  cji  =  0.30,  C22  =  0-29,  and  c \2  =  0.31.  The  corresponding  global  model  constants  are 
Cs  =  0.1,  m  =  8.0,  and  n  =  0.14.  Distribution  of  the  cj2  coefficient,  plotted  in  figure  9,  shows 
pockets  of  nearly  perfect  correlation  close  to  the  cavity  walls  and  near  the  region  of  the  down¬ 
stream  free  shear  layer.  Notice  that  the  poorest  correlations  occur  primarily  where  the  model 
TEV  levels  are  expected  to  be  low.  As  an  example,  see  figure  10a  where  the  highest  TEV  levels 
found  in  the  individual  recirculation  planes  are  shown  for  the  DNS  results  at  completion  of  the 
simulation.  This  figure  also  reveals  that  the  model  predicts  TEV  magnitudes  and  distributions  in 
accordance  with  the  experimental  observations.  As  a  final  note,  the  model  quality  was  also 
checked  against  DNS  results  from  a  higher  Reynolds  number  test  case  (Re  =  7500)  where  the 
flow  is  locally  transitional  [13].  The  highest  TEV  levels  computed  on  the  individual  recircula¬ 
tion  planes  are  shown  in  figure  10b.  Indeed,  the  peak  levels  of  TEV  are  concentrated  within  the 
downstream  free  shear  layer. 


LES  RESULTS 

The  uniform  grid  selected  for  the  LES  computations  at  Re  =  10,000  was  101  x  101  x  81 . 
Based  on  the  turbulent  scales  estimated  by  Koseff  and  Street  [12]  for  the  3D  cavity  flow  at 
Re  =  10,000,  this  grid  provides  a  higher  spatial  resolution  than  that  needed  to  resolve  Taylor’s 
microscale.  Shown  in  figures  11a  and  1  lb  are  snapshots  of  the  velocity  field  at  sample  time 
t  =  6.5.  The  spanwise  velocity  vectors  represent  flow  at  the  x  =  0.77  plane,  while  the  DSE  region 
is  shown  at  various  planes  between  the  TGL  vortices.  Besides  the  comer  vortex,  the  spanwise 
velocity  vectors  show  five  additional  vortex  structures  that  appear  distorted  when  compared  with 
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Figure  7.  Time  Traces  (a)  and  (c)  and  Power  Spectra  (b)  and  (d)  of  Horizontal  and  Vertical  Velocity 
Components  Near  the  Downstream  Free  Shear  Layer  for  Test  Case  Re  =  5000 
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Figure  8.  Quantitative  Comparisons  of  the  DNS  Results  and  Experimental 
Data  [26]  in  Terms  of  the  Centerline  Mean  Velocities  at  the  Cavity  Mid-Span  Plane 
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Figure  9.  Distribution  of  Correlation  Coefficient  c12  Using  Cs  =  0.1 
and  Parameters  m  -  8.0  and  n  =  0.14  in  Turbulence  Length  Scale  Definition. 
Correlation  Computed  from  DNS  Data  Taken  at  Completion  of  Test  Case  Re  =  5000 
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(a)  (b) 

Figure  10.  Turbulent  Eddy  Viscosity  Levels  (Normalized  by  the  Kinematic 
Viscosity)  Computed  from  DNS  Taken  at  Completion  of  Simulation; 

(a)  Re  =  5000  and  (b)  Re  =  7500 


(a)  (b) 


Figure  11.  Velocity  Vectors  Illustrating  (a)  Regions  of  the  Downstream 
Secondary  Eddy  and  (b)  TGL  and  Corner  Vortices  for  LES  Test  Case 
Re  =  10,000  at  Relative  Time  t  =  6.5 
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the  DNS  results  (Re  =  5000),  Breakdown  of  the  TGL  vortex  structure  is  due  to  the  onset  of 
turbulence  within  the  adjacent  downstream  free  shear  layers,  Prasad  and  Koseff  [26]  (SAR  = 
1.0)  and  Koseff  and  Street  [13]  (SAR  =  3.0)  also  reported  a  loss  of  TGL  vortex  structure  at  this 
Re,  but  they  were  unable  to  visualize  the  vortex  flow  patterns  due  to  the  rapid  lateral  dispersion 
of  the  dye  streaks.  The  existence  of  the  vortex  structures  and  their  streamwise  extent  is  further 
illustrated  in  figures  12  and  13  in  the  form  of  stream  function  contours.  Only  results  in  the  lower 
half  of  the  cavity  are  shown.  (The  gridlines  in  the  figures  and  in  the  velocity  vectors  (figure  1  lb) 
help  quantify  the  intensity,  spanwise  size,  and  streamwise  extent  of  each  vortex.)  The  instability 
mechanism  of  TGL  vortex  generation  still  exists  at  this  Re  by  evidence  of  the  vortex  structure 
created  above  the  DSE  between  gridlines  k  =  49  and  k  =  62.  Like  the  DNS  results,  neither  the 
TGL  vortices  nor  the  comer  vortex  interact  with  the  USE.  Thus,  the  DSE  structure  remains 
intact  throughout  the  cavity  span.  On  the  other  hand,  the  irregular  development  of  the  DSE  in 
the  spanwise  direction  suggests  complicated  interactions  between  the  unsteady  effects  of  the 
primary  recirculation  vortex  and  the  TGL  vortices.  As  an  example,  Koseff  and  Street  [17]  and 
Prasad  et  al.  [10]  reported  appearances  of  spiraling  spanwise  motions  within  the  DSE,  which 
they  attributed  to  its  interaction  with  the  local  vortices.  Verification  of  their  observations  and  the 
uncovering  of  other  interactions  are  addressed  next. 


Figure  12.  Stream  Function  Contours  Showing  the  Downstream  and 
Upstream  Secondary  Eddies  of  the  Shear-Driven  Cavity  at 
Re  =  10,000;  Time  t  =  6.5 
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Figure  13.  Stream  Function  Contours  Showing  the  TGL  and  Corner  Vortices  of 
the  Shear-Driven  Cavity  at  Re  =  10,000;  Time  t  =  6.5 


Between  the  spanwise  vortices  shown  in  figure  13,  the  streamlines  form  four  surfaces  that 
give  a  cave-like  impression.  Beneath  these  surfaces,  the  basic  2D  structure  of  the  DSE  develops 
due  to  separation  of  the  primary  recirculation  vortex  from  the  downstream  wall  (see  figure  11a). 
The  vortices  adjacent  to  the  DSE  strongly  influence  its  spanwise  characteristics.  An  attempt  to 
understand  these  complicated  characteristics,  as  well  as  the  other  flow  features  of  the  3D  cavity 
at  Re  =  10,000,  is  illustrated  in  figure  14  at  sample  time  t  =  6.5.  In  this  figure,  five  sets  of  par¬ 
ticle  traces  were  initiated  either  within  the  DSE  region  or  the  outer  extremes  of  the  primary 
recirculation  vortex,  in  particular,  the  first  computational  point  off  the  downstream  wall 
(x  =  0.99).  Also,  each  set  originated  halfway  up  the  downstream  wall  (y  =  0.5)  except  set  5, 
which  started  at  y  =  0.05.  Set  1  contains  three  particle  traces  that  were  initiated  directly  over  the 
center  of  the  large  TGL  vortex  (gridline  number  k  =  10).  After  release,  all  three  particles  were 
entrained  by  the  downwash  flow  of  the  large  vortex.  Their  spiraling  path  then  transversed 
streamwise  only  a  short  distance  before  being  entrained  by  the  upwash  region  of  the  vortex  as 
given  by  the  mid-span  recirculation  flow.  Particles  from  sets  3  and  4  had  a  similar  fate.  Set  3 
initiated  above  and  to  the  left  of  the  TGL  vortex  at  gridline  k  =  28,  while  set  4  originated  be¬ 
tween  the  TGL  vortices  marked  gridlines  k  =  28  and  k  =  36.  Initially,  each  particle  was  con- 
vected  by  the  primary  recirculation  vortex,  which  traced  a  path  that  coincides  with  the  streamline 
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Figure  14.  Five  Sets  of  Particle  Traces  Illustrating  the  Interaction  Among 
the  Primary  Recirculation  Vortex ,  Downstream  Secondary  Eddy,  and  TGL 
and  Corner  Vortices  of  the  LES  Computation  at  Relative  Time  t  =  6.5 

surface  shown  in  figure  12.  As  the  particles  neared  the  cavity  bottom,  they  were  entrained  by  the 
respective  vortex  instead  of  the  DSE.  The  reason  for  this  is  illustrated  by  particle  sets  2  and  5. 
Set  2  started  above  and  to  the  right  of  the  TGL  vortex  at  gridline  k  =  28,  while  set  5  was  centered 
just  above  the  tiny  flow  structure  that  is  positioned  between  gridline  k  =  62  and  the  end  wall.  As 
the  particles  of  set  2  approached  the  cavity  bottom,  they  were  entrained  by  the  DSE.  However, 
because  the  large  adjacent  vortex  induced  a  dominate  spanwise  velocity  component  on  the  DSE, 
the  particles  traced  a  broad  spanwise  spiral  that  turned  quickly  streamwise  once  fully  entrained 
by  the  vortex.  This  spanwise  spiral  within  the  DSE  region  is  also  traced  by  particle  set  5.  From 
these  observations,  we  can  conclude  that  the  streamwise  extent  of  the  TGL  vortex  structures 
shown  in  figure  13  are  sustained  through  two  patterns  of  fluid  entrainment.  Close  to  the  down¬ 
stream  wall,  these  vortices  entrain  fluid  from  the  adjacent  DSE  regions,  which  in  turn  extract 
fluid  from  the  primary  recirculation  vortex.  As  mentioned  earlier,  this  flow  pattern  was  also 
observed  locally  in  the  3D  cavity  in  the  flow  visualization  experiments.  Upstream  from  the  DSE 
region,  however,  the  TGL  vortices  entrain  fluid  from  the  primary  vortex  directly.  Loss  of  the 
vortex  structure  occurs  upstream  once  the  primary  vortex  separates  from  the  cavity  bottom. 
Having  this  understanding,  the  tiny  flow  structure  within  the  DSE  region  between  gridline  k  =  62 
and  the  span’s  end  wall  is  another  vortex  pair.  Creation  of  this  secondary  vortex  pair  is  due  to 
the  opposing  spanwise  viscous  interactions  of  the  comer  vortex,  the  adjacent  TGL  vortex,  and 
the  no-slip  condition  along  the  cavity  bottom.  Above  its  center  position  lies  an  imaginary  sur¬ 
face  within  the  DSE  region  that  demarcates  fluid  entrained  by  the  corner  vortex  apart  from  that 
extracted  by  the  adjacent  TGL  vortex.  The  existence  of  this  flow  surface  is  illustrated  by  the 
right  and  left  particle  traces  in  set  5.  As  shown  in  figure  14,  the  spanwise  spiral  trace  of  each 
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particle  is  in  opposite  directions.  Thus,  these  particles  started  on  opposite  sides  of  the  flow 
surface.  Because  the  adjacent  TGL  vortices  change  size  and  meander  along  the  cavity  floor  over 
time,  the  flow  surfaces  and  corresponding  secondary  vortices  are  not  stationary.  As  a  matter  of 
fact,  flow  visualizations  of  the  numerical  results  showed  extremely  complicated  dynamics  that 
governed  these  tiny  structures.  However,  their  impact  on  momentum  and  energy  transfer 
throughout  the  cavity  was  never  significant. 

Now  the  effects  of  the  symmetry  plane  assumption  on  these  LES  results  will  be  discussed. 
Koseff  and  Street  [12]  noted  that  the  mean  recirculation  flow  in  the  3D  cavity  with  SAR  =  3.0 
became  self-similar  when  the  Reynolds  number  was  increased  to  10,000.  By  examining  the 
comparisons  between  the  computed  centerline  mean  velocity  profiles  and  the  experimental 
results  in  figures  15a  and  15b,  one  can  see  that  this  is  indeed  the  case.  While  the  experimental 
data  are  at  the  mid-span  plane,  the  computed  profiles  are  shown  along  the  centerline  of  the 
recirculation  plane  z  =  0.28  (K15  plane);  zmax  =  1.5.  In  both  profiles,  the  averaged  error  (as 
compared  with  the  data)  is  less  than  2  percent.  Thus,  the  symmetry  plane  assumption  has  a 
negligible  effect  on  the  mean  recirculation  at  the  K15  plane  and  throughout  the  remainder  of  the 
cavity  span.  Unfortunately,  no  centerline  Reynolds  stresses  or  turbulent  velocity  fluctuations  of 
the  3D  cavity  flow  with  SAR  =  3.0  were  published.  However,  data  were  reported  by  Prasad  and 
Koseff  [26]  at  the  mid-span  plane  for  SAR  =  1.0.  Quantitative  comparisons  of  those  data 
and  the  LES  results  at  the  mid-span,  K15,  and  end  planes  are  shown  in  figures  16a  and  16b  in 
terms  of  the  centerline  root-mean-square  (RMS)  velocity  fluctuations.  The  end  plane  profile  in 
each  figure  was  computed  at  recirculation  plane  z  =  1.25.  Both  figures  clearly  show  large  dis¬ 
crepancies  between  the  computed  and  experimental  mid-span  plane  RMS  results.  Conversely, 
RMS  profiles  at  the  K15  and  end  planes  agree  reasonably  well  with  each  other,  as  well  as  with 
the  experimental  mid-span  plane  data.  Furthermore,  the  computed  profiles  extend  the  overall 
observed  trend  that  turbulent  kinetic  energy  is  lost  near  the  cavity  walls  when  the  cavity  SAR  is 
reduced  [17,26],  Hence,  for  this  test  case,  these  observations  along  with  the  experimental  agree¬ 
ment  of  the  computed  mean  recirculation  flow  illustrated  a  localized  effect  of  the  spanwise 
symmetry  plane  assumption  on  the  LES  results. 

As  final  issues  to  address,  do  these  LES  computations  support  the  explicit  understanding 
that  turbulence  dominates  the  flow  characteristics  near  the  downstream  free  shear  layer,  and  at 
this  Reynolds  number  what  are  the  spatial  distributions  of  the  velocity  fluctuations  elsewhere  in 
the  cavity?  Figures  17a  and  17b  show  the  energy  spectrum  of  the  velocity  fluctuations  calcu¬ 
lated  from  a  90-second  sample  record  taken  in  the  vicinity  of  the  downstream  free  shear  layer. 
Both  spectral  profiles  display  an  inertial  subrange  with  about  an  order  of  magnitude  larger 
energy  content  within  the  vertical  fluctuations.  This  result,  as  well  as  the  bandwidth  over  which 
the  inertial  subrange  occurs,  agrees  with  similar  analyses  conducted  on  the  experimental  data 
[12].  In  figures  18a  and  18b,  distributions  of  the  vertical  RMS  fluctuations  and  shear  stress 
component  u'v'  at  the  K15  plane  indeed  disclose  highest  levels  within  the  downstream  free  shear 
layer  region.  However,  one  should  note  that  the  RMS  levels  of  the  vertical  velocity  fluctuation 
are  also  significant  along  the  downstream  wall. 
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Figure  15.  Comparison  of  LES  and  Experimental  Results  [  12]  in  Terms  of 
Centerline  Mean  Velocity  Profiles;  Simulation  Profile  at  K15  Plane  and 
Experimental  Data  Measured  at  Mid-Span  Plane 


Figure  16.  Profiles  of  Centerline  Horizontal  and  Vertical  RMS  Velocity 
Fluctuations;  LES  Profiles  are  at  the  K15,  End,  and  Mid-Span  Planes  While 
Experimental  Data  Are  at  the  Mid-Span  Plane  of  Cavity  with  SAR  =  1:1 
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FREQUENCY 


FREQUENCY 


(a)  (b) 

Figure  17.  Power  Spectra  of  Velocity  Components  Near  the  Downstream 
Free  Shear  Layer;  (a)  Horizontal,  (b)  Vertical 


Figure  18.  Distribution  of  the  (a)  Vertical  RMS  Fluctuations  and 
(b)  u'v'  Shear  Stress  Component  at  the  K15  Plane 
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CONCLUSIONS 


Direct  numerical  simulation  (DNS)  and  large-eddy  simulation  (LES)  methodologies  are 
used  to  study  the  unsteady  and  turbulent  characteristics  of  the  three-dimensional  shear-driven 
cavity  flow  at  Reynolds  numbers  5000  and  10,000.  Based  on  the  numerical  results,  the  follow¬ 
ing  conclusions  about  this  flow  are  offered: 

1.  At  Re  =  5000  (DNS  results),  the  flow  is  laminar,  but  the  three-dimensionality  and 
unsteadiness  severely  disturbs  the  basic  structure  of  the  classic  recirculation  flow  features.  The 
TGL  vortices  change  rapidly  in  size  and  they  meander  only  locally.  We  predicted  nine  TGL 
vortex  pairs  spanning  the  cavity  bottom,  with  one  straddling  the  mid-span  plane.  The  flow 
remains  fully  deterministic.  Since  the  onset  of  turbulence  occurs  between  Reynolds  numbers  of 
6000  and  8000  [12],  these  results  provided  a  reasonable  means  for  determining  the  global  param¬ 
eters  in  the  SGS  length  scale  definition. 

2.  At  Re  =  10,000  (LES  results)  the  instability  mechanism  for  TGL  vortex  generation 
still  exists,  but  the  vortices  themselves  have  now  become  distorted  due  to  the  onset  of  turbulence 
within  the  downstream  free  shear  layer.  Their  physical  characteristics  throughout  the  3D  cavity 
vary  randomly.  The  combined  effects  of  the  primary  recirculation  vortex  and  the  TGL  and 
corner  vortices  cause  a  complicated  irregular  development  of  the  DSE.  Near  the  downstream 
wall  the  TGL  vortices  extract  fluid  from  the  downstream  eddy  region,  which  concurrently  en¬ 
trains  fluid  from  the  primary  recirculation  vortex.  Upstream  of  the  downstream  eddy  region,  the 
TGL  vortices  entrain  fluid  directly  from  the  primary  vortex.  Lastly,  secondary  vortex  pairs  are 
created  intermittently  within  the  DSE  region  due  to  the  viscous  interactions  among  the  adjacent 
larger  vortices  and  the  cavity  floor. 
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